home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / c / jpl_c.zip / INVERF.C < prev    next >
Text File  |  1986-05-18  |  6KB  |  197 lines
  1. /* 1.1  10-08-85                        (inverf.c)
  2.  ************************************************************************
  3.  *            Robert C. Tausworthe                *
  4.  *            Jet Propulsion Laboratory            *
  5.  *            Pasadena, CA 91009        1985        *
  6.  ************************************************************************/
  7.  
  8. #include "defs.h"
  9. #include "stdtyp.h"
  10. #include "errno.h"
  11. #include "mathcons.h"
  12.  
  13. /*----------------------------------------------------------------------*
  14.  
  15.     This routine uses the approximations given in J. M. Blair,
  16. et. al., "Rational Chebyshev Approximations for the Inverse Error
  17. Function," Math Comp. Vol. 30, No. 136, Oct. 1976, pp. 827-830, with
  18. the appendices given in separate microfiche.
  19.  
  20.  *----------------------------------------------------------------------*/
  21.  
  22. LOCAL int dP[4] = {7, 8, 11, 8};    /* degrees of P(x)        */
  23. LOCAL int dQ[4] = {7, 8, 10, 9};    /* degrees of Q(x)        */
  24.  
  25. LOCAL double P[4][12] =        /* numerator polynomials        */
  26. {
  27. /* P[0][0] */    -.30866886527764497339e2, /* Table 19 of Ref. below    */
  28. /* P[0][1] */    0.20652786402942339589e3, /*    |x| <= 0.75        */
  29. /* P[0][2] */    -.52856770835093823310e3,
  30. /* P[0][3] */    0.64880509544036249088e3,
  31. /* P[0][4] */    -.39205509901971391289e3,
  32. /* P[0][5] */    0.10706278097770074402e3,
  33. /* P[0][6] */    -.10303488455439678272e2,
  34. /* P[0][7] */    0.15641510821923860000e0,
  35. /* P[0][8] */    0.0,
  36. /* P[0][9] */    0.0,
  37. /* P[0][10] */    0.0,
  38. /* P[0][11] */    0.0,
  39.  
  40. /* P[1][0] */    0.94897362808681080020e-2, /* Table 39 of Ref. below    */
  41. /* P[1][1] */    -.24758242362823355485e0,  /* 0.75 <= x <= 0.9375    */
  42. /* P[1][2] */    0.25349389220714893916e1,
  43. /* P[1][3] */    -.12954198980646771502e2,
  44. /* P[1][4] */    0.34810057749357500873e2,
  45. /* P[1][5] */    -.47644367129787181802e2,
  46. /* P[1][6] */    0.29631331505876308122e2,
  47. /* P[1][7] */    -.64200071507209448654e1,
  48. /* P[1][8] */    0.21489185007307062000e0,
  49. /* P[1][9] */    0.0,
  50. /* P[1][10] */    0.0,
  51. /* P[1][11] */    0.0,
  52.  
  53. /* P[2][0] */    0.19953288964537210824e-5, /* Table 60 of Ref. below    */
  54. /* P[2][1] */    0.21273702631785953343e-3, /* 0.9375 <= x <= 1 - 10^(-100) */
  55. /* P[2][2] */    0.58975595952407247651e-2,
  56. /* P[2][3] */    0.59481901452735809123e-1,
  57. /* P[2][4] */    0.31328289030932667506e0,
  58. /* P[2][5] */    0.13630199956442260990e1,
  59. /* P[2][6] */    0.34152815205652930673e1,
  60. /* P[2][7] */    0.30184181468933606839e1,
  61. /* P[2][8] */    0.20842433546207339433e1,
  62. /* P[2][9] */    0.85545635026743499993e0,
  63. /* P[2][10] */    0.40273918408712893132e-2,
  64. /* P[2][11] */    -.15196139115744716810e-3,
  65.  
  66. /* P[3][0] */    .45344689563209398449e-11, /* Table 82 of Ref. below    */
  67. /* P[3][1] */    .46156006321345332510e-8,  /* 1 - 10^(-100) <= x    */
  68. /* P[3][2] */    .12964481560643197452e-5,  /*     <= 1 - 10^(-10000)    */
  69. /* P[3][3] */    .13714329569665128933e-3,
  70. /* P[3][4] */    .60537914739162189689e-2,
  71. /* P[3][5] */    .11279046353630280004e0,
  72. /* P[3][6] */    .82810030904462690215e0,
  73. /* P[3][7] */    .19507620287580568829e1,
  74. /* P[3][8] */    .69952990607058154857e0,
  75. /* P[3][9] */    0.0,
  76. /* P[3][10] */    0.0,
  77. /* P[3][11] */    0.0
  78. };
  79.  
  80. LOCAL double Q[4][11] =        /* denominator polynomials        */
  81. {
  82. /* Q[0][0] */    -.28460290173882339383e2,    /* Table 19        */
  83. /* Q[0][1] */    0.20518924149238530630e3,
  84. /* Q[0][2] */    -.57617125090127638064e3,
  85. /* Q[0][3] */    0.79669388170563770334e3,
  86. /* Q[0][4] */    -.56509305564023424022e3,
  87. /* Q[0][5] */    0.19450320483954087700e3,
  88. /* Q[0][6] */    -.27371852306002662877e2,
  89. /* Q[0][7] */    1.0,
  90. /* Q[0][8] */    0.0,
  91. /* Q[0][9] */    0.0,
  92. /* Q[0][10] */    0.0,
  93.  
  94. /* Q[1][0] */    0.67544512778850945940e-2,    /* Table 39        */
  95. /* Q[1][1] */    -.18611650627372178511e0,
  96. /* Q[1][2] */    0.20369295047216351160e1,
  97. /* Q[1][3] */    -.11315360624238054876e2,
  98. /* Q[1][4] */    0.33880176779595142684e2,
  99. /* Q[1][5] */    -.53715373448862143348e2,
  100. /* Q[1][6] */    0.41409991778428888715e2,
  101. /* Q[1][7] */    -.12831383833953226499e2,
  102. /* Q[1][8] */    1.0,
  103. /* Q[1][9] */    0.0,
  104. /* Q[1][10] */    0.0,
  105.  
  106. /* Q[2][0] */    .19953210379374212953e-5,    /* Table 60        */
  107. /* Q[2][1] */    .21274156963404084598e-3,
  108. /* Q[2][2] */    .59037062023731354671e-2,
  109. /* Q[2][3] */    .59959150110861092334e-1,
  110. /* Q[2][4] */    .32318083080817836442e0,
  111. /* Q[2][5] */    .14378337804749344527e1,
  112. /* Q[2][6] */    .37644571508257969664e1,
  113. /* Q[2][7] */    .44081435698143841173e1,
  114. /* Q[2][8] */    .42508710497182804606e1,
  115. /* Q[2][9] */    .22127469427969785343e1,
  116. /* Q[2][10] */    1.0,
  117.  
  118. /* Q[3][0] */    .45344687377088206782e-11,    /* Table 82        */
  119. /* Q[3][1] */    .46156017600933592558e-8,
  120. /* Q[3][2] */    .12964671850944981712e-5,
  121. /* Q[3][3] */    .13715891988350205065e-3,
  122. /* Q[3][4] */    .60574830550097140404e-2,
  123. /* Q[3][5] */    .11311889334355782064e0,
  124. /* Q[3][6] */    .84001814918178042918e0,
  125. /* Q[3][7] */    .21238242087454993541e1,
  126. /* Q[3][8] */    .15771922386662040545e1,
  127. /* Q[3][9] */    1.0,
  128. /* Q[3][10] */    0.0
  129. };
  130.  
  131. /************************************************************************/
  132.     double
  133. inverf(x)        /* inverse of erf(x), which see.        */
  134.  
  135. /*----------------------------------------------------------------------*/
  136. double x;
  137. {
  138.     int interval, sign;
  139.     double y, inverfc(), ratfun();
  140.     LOCAL double bias[2] = {0.5625, 0.87890625};
  141.  
  142.     sign = (x < 0.0 ? 1 : 0);
  143.     if (sign)
  144.         x = -x;
  145.     if (x >= 1.0)
  146.     {    errno = EDOM;
  147.         return (sign ? -INFINITY : INFINITY);
  148.     }
  149.     if (x > 0.9375)
  150.     {    y = inverfc(1.0 - x);
  151.         return (sign ? -y : y);
  152.     }
  153.     if (x < 0.75)
  154.         interval = 0;
  155.     else    interval = 1;
  156.     y = x * ratfun(x*x - bias[interval], P[interval], Q[interval],
  157.         dP[interval], dQ[interval]);
  158.     return (sign ? -y : y);
  159. }
  160.  
  161. #define BOUNDARY    0.065901028    /* 1 / sqrt(100 * log(10))    */
  162.  
  163. /*\p*********************************************************************/
  164.     double
  165. inverfc(x)    /* returns double inverse of erfc(y)            */
  166.  
  167. /*----------------------------------------------------------------------*/
  168. double x;
  169. {
  170.     int sign, interval;
  171.     double y, log(), fabs(), inverf(), sqrt(), ratfun();
  172.  
  173.     y = fabs(1.0 - x);
  174.     if (y >= 1.0)
  175.     {    errno = EDOM;
  176.         return (x > 2 ? -INFINITY : INFINITY);
  177.     }
  178.     if (x > 1.0)
  179.     {    x = 2.0 - x;
  180.         sign = 1;
  181.     } else    sign = 0;
  182.     if (x IS 0.0)
  183.     {    errno = ERANGE;
  184.         return (sign ? -INFINITY : INFINITY);
  185.     }
  186.     if (x > 0.0625)
  187.     {    y = inverf(1.0 - x);
  188.         return (sign ? -y : y);
  189.     }
  190.     y = 1.0 / sqrt(-log(x));
  191.     if (y > BOUNDARY)
  192.         interval = 2;
  193.     else    interval = 3;
  194.     y = ratfun(y, P[interval], Q[interval], dP[interval], dQ[interval])/y;
  195.     return (sign ? -y : y);
  196. }
  197.